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We present an explicit solver of the three-dimensional screened and unscreened Poisson's equation 
which combines accuracy, computational efficiency and versatility. The solver, based on a mixed 
plane-wave / interpolating scaling function representation, can deal with any kind of periodicity 
(along one, two, or three spatial axes) as well as with fully isolated boundary conditions. It can 
seamlessly accommodate a finite screening length, non-orthorhombic lattices and charged systems. 
This approach is particularly advantageous because convergence is attained by simply refining the 
real space grid, namely without any adjustable parameter. At the same time, the numerical method 
features 0(N log N) scaling of the computational cost (N being the number of grid points) very much 
like plane-wave methods. The methodology, validated on model systems, is tailored for leading-edge 
computer simulations of materials (including ab initio electronic structure computations), but it 
might as well be beneficial for other research domains. 



I. INTRODUCTION 

Poisson's equation (screened or not) is involved in a 
large variety of problems in physics and chemistry as 
well as in engineering. Therefore, there is a quite strong 
motivation for developing efficient and accurate solving 
methods. 

As far as electrostatics is concerned, the three- 
dimensional screened Poisson's equation is written as 
follows (in Gaussian units): 

(V 2 -til)V(x,y,z) = -4irp(x,y,z), (1) 

where p(x, y, z) represents a continuous electric charge dis- 
tribution (the input of the problem at hand), V(x, y, z) is 
the electrostatic potential (the output), and /x represents 
the reciprocal screening length as defined, for instance, 
in the Debye-Hiickel or Thomas- Fermi approximations. 
In the special case fiQ = 0, Eq. (1) reduces to the usual 
Poisson's equation. 

Any method aiming at providing a solution to Eq. (1) 
has to deal with boundary conditions (BC), which in 
general can be cither periodic or free (otherwise referred 
to as "isolated" or "open"') along each of the three direc- 
tions x,y,z. In the case of fully periodic BC, the most 
natural (and efficient) approach to the problem is that of 
the reciprocal space treatment. It amounts to expanding 
both the density and the potential as superpositions of 
plane waves (Fourier series), following which Eq. (1) be- 
comes algebraic in the Fourier components of p and V. 
This equation is readily solved and the result is finally 
transformed back into real space. Forward and backward 
transformations are carried out via Fast Fourier Trans- 
form (FFT) , hence the overall computational scaling of 
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the method with respect to the number N of grid points 
is a rather appealing (TV log TV). 

Owing to the above mentioned advantages, many at- 
tempts have been made to tackle also free BC or mixed 
free/periodic BC within a nominally fully periodic frame- 
work. In the simplest implementations, the simulation 
box is artificially enlarged by vacuum padding, such as to 
suppress the spurious Coulombian interaction among peri- 
odic replicas ("super-cell" approximation). Yet, as a result 
of the the long-range nature of the Coulomb interaction, 
there are situations in which the simulation box ought 
to be made unfeasibly large, in particular for charge dis- 
tributions exhibiting significant multipolar terms. More- 
over, a non-zero net charge in the primitive cell would 
yield a divergent total electrostatic energy when infinitely 
replicated along every direction, unless a compensating 
artificial uniformly charged background ("jellium") is in- 
troduced, see e.g. Refs. [I] and [2]. Another option which 
has been put forward consists in cutting off the long tail 
of the Coulomb interaction beyond a spherical region in 
real space, the radius of which is adequately chosen. Cor- 
respondingly, the reciprocal space components of the bare 
Coulomb potential are multiplied by screening functions, 
known analytically for all types of periodicity [3, 4]. In a 
series of papers [5-7] , the screening function formalism 
was combined with the explicit break-up of the short- 
and long-range components of the Coulomb interaction, 
as also done in the context of (smooth particle-mesh) 
Ewald summation techniques, whereas in Refs. [8, 9] the 
errors induced by the periodic images are alleviated by 
introducing a corrective potential. 

Although we acknowledge that much remarkable work 
has been done on the subject, providing a thorough review 
would go beyond the scope of the present study. We 
therefore refer the reader to the original literature and 
proceed by presenting our approach, which differs from 
those mentioned above in that convergence is attained 
with no adjustable parameter and thus it aims at being 
fully generic. 
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This paper builds upon Refs. [10, 11], where a novel 
method for solving the unscreened Poisson's equation 
with free and surface-like BC was first presented. Such a 
method is direct (rather than iterative) in that the solution 
along the isolated directions is found in its integral form 
by using the Green's function method. For instance, in 
the case of a fully isolated system (or "cluster-like"), 

V(r) = Awf dV G( Mo ; \r - ?\)p(?), (2) 

where f= (x,y,z). Homogeneous Dirichlet BC (V = 
at \r\ — > oo) along the isolated directions are explicitly 
enforced by the selection of the Green's function. 

The method has been in use for a few years in a number 
of ab-initio codes, namely ABINIT [12, 13], BigDFT 
[14, 15], CP2K [16] (see also Ref. [17] for a recent ap- 
plication thereof), OCTOPUS [18, 19] and has proven to 
be highly efficient and accurate in every application at- 
tempted to date. It is based on a mixed plane-wave / 
interpolating scaling function (ISF) representation of p 
and V which allows to model any sort of periodicity in 
the most natural, clean and mathematically rigorous way. 
Clearly, periodic (isolated) directions are represented in 
terms of plane waves (interpolating scaling functions). 
ISFs - arising in the wavelet theory [20, 21] - enjoy sev- 
eral properties which make them superior to other basis 
sets. For instance, the representation in terms of m-th 
order ISFs make the first m moments of the continuous 
and discrete charge distributions coincide [11]. As a con- 
sequence the representation is definitely faithful (other 
than handy) , since the different moments of the charge 
distribution capture the major features of the potential. 
Moreover, ISFs arc genuinely localized due to their com- 
pact support (the length of which is equal to 2m) and 
endowed with the so-called "refinement relations" which 
easily allow to switch from a representation on a grid with 
spacing h to a doubly refined grid with spacing h/2. 

We have extended the previous implementation to 
account also for screening, for the case of periodicity 
along only one direction ("wire-like BC") and for non- 
orthorhombic cells, this investigation providing a detailed 
account of such improvements. The inclusion of such 
new functionalities is motivated by the strong theoretical, 
experimental and technological interest in the characteri- 
zation of nanostructured materials (among which polar 
nanorods, see e.g. Ref. [22] and references therein), since 
solving Poisson's equation is only one of the many steps 
involved in state-of-the-art computer simulations and is 



repeated several times. Moreover, in the context of Kohn- 
Sham (KS) density functional theory (DFT) and exten- 
sions thereof, there are quantities which are computed 
via convolution integrals very similar to that in Eq. (2): 
for instance, the exact exchange term arising within those 
generalizations of KS-DFT employing orbital-dependent 
or hybrid functionals (see [23] and references therein), 
or the coupling-matrix in time-dependent DFT [24] . In 
this respect, the electrostatic problem of concern here 
provides the paradigm for many other computations, even 
well beyond the scope of electrostatics. 

Regarding the possibility of accounting also for screen- 
ing, we note that this novel feature might for example 
be used to solve the Schrodinger equation iteratively (see 
e.g. Ref. [25]). In fact, one can exploit the formal anal- 
ogy between Eq. (1) and the Schrodinger equation which 
becomes apparent if the latter is written in the following 
fashion: 

(|V-|£|)|V>> = F|V>>, (3) 

where \?p) represents a bound eigenstate with negative 
energy (E < 0). 

The next sections are structured as follows: we first 
present our solution method for free, wire-like and surface- 
like BC. We then discuss the accuracy of the proposed 
solver by reporting a collection of numerical benchmarks. 
We conclude by highlighting the benefits of using the 
present methodology. 



II. FREE BOUNDARY CONDITIONS 

In the case of free BC, the Green's function which has 
to be plugged into Eq. (2) is 



since 

(^ + ll-d)G{^r) = -^\r), (5) 

where r = \r\. Along the same lines as Ref. [11], both 
the charge distribution and the electrostatic potential are 
expanded in terms of ISFs, here denoted by </>: 
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where h^ x y z y and N^ x y z y represent the (uniform) grid respectively. Since the function (j) is, by construction, 
spacing and the number of grid points along each direction, 
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such that 4>(j) = Sj.o, Vj S Z, the expansion coefficients 
Pjsedydn are readily found to be: 

Pj*,jyj. = P( h xjx, h y j yi h z j z ). (7) 

In other words, the expansion coefficients coincide exactly 
with the values of p(x,y,z) on a uniform grid. In this 
respect, the ISF representation appears genuinely tailored 
for numerical studies, where the whole available infor- 
mation reduces to knowing the values at the grid points. 
Clearly, the grid spacing has to be chosen adequately, 
depending on the typical spatial scales over which the 
density distribution exhibits significant variations. The 
underlying mathematics then assures that the moments 
built upon the discrete charge distribution coincide with 
those of the continuous charge distribution up to order 
m — 1, where m is the order of the ISF (see Ref. [11] for 



the proof): 

ih 3 h kh PW = J d 3 rx h y l -z l3 p(r) (8) 

if < li, I2, 13 < Tn. A representation analogous to that 
in Eq. (6) can also be given for the potential V, where 

V j*,jy,j* = v ( h xjx, hyj y , h z j z ) (9) 

replaces Pj x j j 2 - As a consequence of the chosen repre- 
sentation, the convolution integral in Eq. (2) is more con- 
veniently expressed in Cartesian coordinates, and upon 
plugging Eq. (6) into Eq. (2) we obtain the following 
equation in discrete form: 



N x N y 7V S 

VjbJvJ* = 4>*h x hyh z K(j x - ./'../., - f y ,jz - ./':/'••) Pj'^j^, (10) 

j'4=oj'i=oii=Q 

where 

K (jx,jy,jz;p) = J dudvdwG(p; \J[h x {j x - u)] 2 + [h y (j y - v)] 2 + [h z (j z - w)] 2 ) <j)(u)(f>(v)(f>(w). (11) 
is the convolution kernel. , 



The numerical evaluation of Eq. (11) would be too oner- 
ous if performed directly. The computational effort can be 
drastically reduced by expressing the Green's function as 
a linear combination of Gaussian functions, as the integral 
would become separable along x, y, z and the resulting 
ID integrals can be evaluated very efficiently. 

We proceed by approximating the Green's function as 

ft 

where a^i/J-) and u>k(p) are determined so as to minimize 
the error on a given range of r. More specifically, we 
found out the following approximation, 

— ^]>> fc e-^ 2 , (13) 
x — ' 

ft=i 

with satisfactory accuracy for any x £ [10 _9 ,33] (see 
Fig. 1). The best fit was performed using the Levenberg- 
Marquardt algorithm [26, 27] with 136 Gaussian functions, 
the afe- values ranging between 10 -3 and 10 19 . The actual 
Green's function can therefore be written as in Eq. (13) 
with 

ui k (fj.) = atkifi) = M 2 «fc- (14) 



In the unscreened case (po = 0) we use the Gaussian fit of 
the 1/r function which was already proposed in Ref. [11] 
and that we recall as being affected by an error < 10~ 8 
for any r £ [1(T 9 ,1]. 

The convolution kernel can be written as follows: 

K(j x ,jy,jz;p) = Y UJk ^ II I ( a k(p)hl,ji) , 

(15) 

where 

I(a;j) = J dte- a W<f>(t), (16) 

and can be computed by evaluating (N x + N y + N z )Ng 
ID integrals - Nq being the number of Gaussian functions 
- hence at a much lower cost than Eq. (11), which would 
require the computation of N x N y N z 3D integrals, instead. 

We point out that the numerical evaluation of Eq. (16) 
is performed using the same method described in Ref. 
[11], which exploits the refinement relations fulfilled by 
the ISFs and yields an accuracy as high as the machine 
precision, even for the narrowest Gaussians. To this 
effect, ISFs show their superiority over other basis sets 
(cf. e.g. the explicit method laid out in Ref. [28], where 
a Gaussian approximation similar to ours is carried out 
within a discrete variable representation approach) . 
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Figure 1: Accuracy of the approximation of the function 
e~ x /x with 136 Gaussians used in the solution of the screened 
Poisson's equation for the case of free BC. The range of the 
independent variable x is [10 _9 ,33]. We plot both the absolute 
and the relative error because the latter is a better indicator 
close to the origin (where the fitted function takes on very 
large values), while the former is a reliable signature of the 
goodness of the fit towards the opposite end. Note that at 
x = 33 the function e~ x jx is already smaller than the machine 
precision. 



III. WIRE-LIKE BOUNDARY CONDITIONS 

We now consider a system which is periodic along the z 
direction (with period equal to L z ) and isolated over the 
xy-plane. We can hence expand the continuous charge 
density distribution as a sum over its Fourier components 
along z: 



p(x,y,z) =53' 



-2-Ril 



z p v A x -,y) > 



(17) 



noting that in this case the Fourier coefficients are solely 
functions of x and y. After expanding the electrostatic 
potential in a similar manner, the screened Poisson's equa- 
tion yields the following relation between the potential's 
reciprocal space components and those of the density: 



[d 2 x + d 2 -p 2 -p 2 Pz ]V Pz (x,y) 



-Awp Pz (x,y) , (18) 



where p Pz = 2irp z /L z . The symmetry of the problem 
suggests writing the Green's function of Eq. (18) in cylin- 
drical coordinates: 



dr 2 



ld_ 
r dr 



G(p;r) = -S^(r), (19) 



where r = \fx 2 + y 2 and p 2 
Eq. (19) is given by 



Mo 



The solution of 



G(p;r) 



JLfKoOur) M>0 
2ir I - log(r) p = ' 



(20) 



where Ko is the zero-th order modified Bessel function 
of the second kind. We then express the 2D Fourier 
components of both density and potential in terms of ISFs, 
thereby completing the required steps towards the mixed 
plane-wave / ISF representation for the case investigated 
here: 



jx=0 jy=0 



W/, h )"{h„ jy 1 



( 21 ) 

where h x and h v are the grid spacings along the non- 
periodic directions. Combining Eq. (18) with Eq. (17) 
and Eq. (21) one obtains 



V u,WP, = 4nh x hy x 



(22) 



where the kernel K(j x ,j y ; p) is very similar to that in Eq. 
(11), except that in (22) the integral is restricted to the 
non-periodic directions. Within the approximation (12) 
the kernel elements can be evaluated as in Eq. (15) with 
i e {x,y}. 

Clearly, the accuracy of the whole method depends on 
the accuracy of the Gaussian approximation of the Green's 
function. Note that close to the origin the function Ko(x) 
behaves like log(a;), whereas it decreases exponentially for 
large arguments: 



Ko(a:) 




log(x) x — > 
x — > oo 



(23) 



We were able to approximate Kq(x) as a sum of Gaussians 
(we denote the fitting parameters by {ui' k , d' k }, cf. Eq. (13)) 
with an accuracy better than 10~ 10 in the range [10 -9 , 30], 
where the upper bound for x was chosen by realizing that 
the value attained by Kq(x) at x — 30 is ~ 2.1 x 10 -14 , 
i.e. already comparable with the machine precision. The 
function fitting was carried out using the same algorithm 
as before, but with 144 Gaussians, with a' k - values ranging 
between 10~ 7 and 10 20 . We used the very same 144 
Gaussians (a' k ' = a' k ) but different relative weights (ui k ) 
also to approximate log(a;), achieving an accuracy of 10 -8 
in the range [10 -9 , 1] (see Fig. 2). 

Therefore, for any p > the convolution kernel is the 
same as in Eq. (15) with i G {x,y}, u>k{p) = w^. and 
otk{p) = Oi' k p 2 , whereas for p — we exploit the scaling 
properties of the logarithm to adapt the best fit obtained 
for x £ [0, 1] to any r 6 [0, L). We write it explicitly for 
sake of clarity: 



K(j x ,j y ;p = 0) = log(L) + 



(24) 
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k h 2 - 7 



where L = sJ[h x {N x + 2m)] 2 + [h y (N y + 2m)} 2 is delib- 
erately chosen so that r/L < 1. In fact, h{. 
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where i,j € {x, y} (i.e. the periodic directions with period 



equal to L x 
tensor, 



L y ) and g lJ is the 2D contra- variant metric 



sin 2 a 



1 

- cos a 



cos a 
1 



(27) 



a being the angle between the x and y unit vectors. 

We point out that in the case of surface-like BC there 
is no Gaussian approximation of the Green's function 
involved in the procedure. Consequently, in this case the 
accuracy of the method is limited only by the machine 
precision. 



NUMERICAL RESULTS 



Figure 2: Accuracy of the approximation of the Green's 
function with 144 Gaussians as used in the solution of the 
screened Poisson's equation for the case of wire-like BC. The 
range of the independent variable x is [10 -9 , 30] for the function 
K (x) and [1(T 9 ,1] for log(x). 



the box size along {x, y}, i.e. the maximum range covered 
by the charge distribution along each axis; 2m hi x , y y is 
the extent of the ISF definition domain. The convolution 
of any charge distribution function with an ISF can be 
non-zero at most over a rectangular domain of diagonal 
length L. 

In order to find the solution for the electrostatic poten- 
tial in real space, we first compute the Fourier coefficients 
of the density (pj x ,j v ; P: ,) through a ID FFT along the 
periodic direction z. The corresponding quantities for 
the potential are then obtained by calculating the con- 
volutions in Eq. (22) via a zero-padded FFT procedure 
[29]. Finally, the potential is transformed back into real 
space along the z direction. Let us notice that real-to- 
complex FFTs can be used instead of complex-to-complex 
FFTs, since all the quantities are real and the kernel is 
symmetric. 



IV. SURFACE-LIKE BOUNDARY CONDITIONS 

While referring the reader to Ref. [10] for a more ex- 
tensive description of our treatment of surface-like BC, 
in the following we just discuss our methodological im- 
provements . In particular, we are now able to allow also 
for the screening and for monoclinic lattices, simply by 
redefining the p? in the equation relating the 2D Fourier 
components of density and potential, 



dz 2 



P 2 I V, 



Px,Py 



(z) = -4tt P Px , Ph (z), 



as follows: 



Mo 



4tt 2 5>^ 

hi 



Li L i 



(25) 



(26) 



In order to measure the accuracy of our method, we used 
several test charge distributions for which the Poisson's 
equation is exactly solvable, and compared the approx- 
imate numerical solution against the exact one. In the 
following the accuracy is given in terms of the infinity 
norm: 



lerrHoo = max 

3& :Jy >3z 



3x i3y i3 z 



v. 



(exact) I 



3x,3y,3z 1 



(28) 



Without loss of generality, all tests were run on a cubic 
simulation box (L x = L y = L z = L), with a uniform grid 
and an equal number of points along each direction. 

In particular, in the case of free BC we used a Gaussian 
density distribution, p(r) — A exp [— r 2 /(2cr 2 )] . As our 
procedure relies on a convolution - see Eq. (2) - from which 
no divergence can arise (as long as p(r) is regular), we 
expect that the numerical solution is regular everywhere. 
The latter statement offers an unambiguous prescription 
for fixing the integration constants of the analytic exact 
solution. The latter is eventually found to be as follows: 



V(r) = A(V2ttct) 
r 



2r 



(29) 



l , [( :|--^ + ^)-^erfcf-L + ^ 
V2a V2 J \V2a V2 



The results of our tests - with A such that V(r — > 0) = 
1 - are reported in Fig. 3 for all the ISF supported in 
our code but no screening and in Fig. 4 for 16th order 
ISF and selected values of /i 2 , covering four orders of 
magnitude. We point out that, owing to the chosen value 
of cr, the value of the density on the simulation box faces 
is pir = L) = 2.21 x 10~ 12 bohr -3 , hence smaller than 
our solver's accuracy. In order to study the influence 
of the box size on accuracy, we performed several runs 
with different box sizes, computed the Hartree energy, 
and compared the result to the exact Hartree energy 
corresponding to an infinite box size and unitary monopole 
(q = j d 3 rp(r) = 1): 



E 



(exact) 



1 



d 3 fp(r)V(r) 



1 

2^V 



(30) 
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Figure 3: Accuracy test for the case of free/isolated boundary 
conditions in the absence of screening (m is the order of the 
ISF, h the grid spacing). 



Figure 4: Accuracy test for the case of free/isolated boundary 
conditions in the presence of screening (m = 16 is the order 
of the ISF, h the grid spacing). 



Free Boundary Conditions - p. = 0, m = 16 



Throughout this second set of tests the value of A was 
chosen so as to keep q — 1. The results are reported 
in Fig. 5, together with the value of the charge density 
distribution at the box faces. We can observe that there 
is absolutely no need to enlarge the box size beyond the 
reference size, as the latter is already large enough to 
capture all the features within the reach of our solver. 
On the other hand, accuracy decreases on reducing the 
box size, as the ideal free BC p(r — > L) — > becomes 
increasingly violated. Nevertheless, the attained accuracy 
remains interesting over a broad range of box sizes even 
smaller than the reference one, at variance with plane- 
wave methods which would fail our test, firstly because of 
the presence of a non-zero monopole and secondly because 
of the aliasing due to insufficiently large box. 

We remark that, as opposed to what claimed in Ref. 
[30] in relation to Refs. [10, 11] (where the method de- 
ployed here was first proposed), our solver is actually 
reliable without any resort to the so-called "minimum 
image convention", namely without rendering the input 
density charge distribution periodic within the simulation 
box. 

In the case of surface-like and wire-like BC, we consider 
a charge density distribution obtained by applying explic- 
itly the screened Poisson's operator on an exact potential 
written as 

V{x,y,z) = -4* H f {P , I} (i;Li) , (31) 

iG{x,y,z} 

where each of the ad hoc functions /'s entering the product 
is either periodic, 



fp(x;L) = exp 



cos 2ir 



-L/2 
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(32) 
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Figure 5: Influence of the (cubic) simulation box size on the 
accuracy of the Hartree energy. L stands for the box size, 
whereas £ re f stands for the box size for which p(r = L) = 



2.21 x 1CT 
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or localized, 

fi(x;L) 
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(33) 



depending on the intended BC. Results are shown in Figs. 
6, 7, 8, 9 and indicate a very good overall convergence 
rate. 

In the case of wire- like BC, we furthered our tests by 
choosing a 2D Gaussian charge distribution, 



p(r, z) = p(r) 



-kr A 



(34) 



where the charge density along z is implicitly set to unity. 
The corresponding exact potential to be used as reference 
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Figure 6: Accuracy test for the case of surface- like boundary 
conditions in the absence of screening (m is the order of the 
ISF, h the grid spacing). 



Figure 8: Accuracy test for the case of wire-like boundary 
conditions in the absence of screening (m is the order of the 
ISF, h the grid spacing). 
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Figure 7: Accuracy test for the case of surface-like boundary 
conditions in the presence of screening (m = 16 is the order 
of the ISF, h the grid spacing). 



is 



V(r) = 4tt 



Ei(— kr 2 ) — log(r 2 ) 
4k 



(35) 



where Ei(a;) is the exponential integral function. On 
deriving Eq. (35), integration constants were fixed unam- 
biguously by following the same criteria which led to Eq. 
(29). In particular, one integration constant guarantees 
the regularity of the potential at the origin, 



lim V(r) = j[ie 

r-s-0 k 



•log(fc)] 



(36) 



where 7b is the Euler-Mascheroni constant, while the 
second integration constant (an additive constant) is set 
to zero so that V(r) ~ — log(r) as r — > oo, in accordance 
with the behavior of the Green's function. This test case 
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Figure 9: Accuracy test for the case of wire-like boundary 
conditions in the presence of screening (m = 16 is the order 
of the ISF, h the grid spacing). 



is relevant to our studies for a two-fold reason. Firstly, 
it is typically out-of-reach for plane-wave methods, since 
the density distribution in Eq. (34) exhibits a non-zero 
monopole. Secondly, it allows to probe the goodness of the 
Gaussian approximation of the log (2;) function involved 
in Eq. (20) within the approximation (12). The charge 
distribution being constant along z, the only non-trivial 
term in the Fourier expansion along z is the zero-mode 
(fi Pz — 0) and, in case the screening is absent (/xq = 0), 
only the log-branch of the Green's function (20) plays a 
role in the computation. We have also analyzed the rate 
of convergence towards the exact Hartree linear energy 
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Figure 10: Accuracy test for the case of wire- like boundary 
conditions (WBC) with monopolar charge density distribution 
(m is the order of the ISF, h the grid spacing). 



density, 



e (exact) ^ 1 I d 2 Fp ( FtZ)V ft z)= (37) 



tt 2 / k 

W*V E + log 2 



(38) 



so as to have a further confirmation that the great ac- 
curacy in the numerically evaluated electrostatic poten- 
tial ensures the reliable computation of derived physical 
quantities. The results are shown in Figs. 10-11 and are 
definitely good. We also include in Fig. 12 a plot of the 
charge density distribution (magnified by a factor 10) 
and the corresponding potentials obtained at different 
/io's. We observe that for /iq = the potential does not 
fall to zero for increasing r, whereas it tends to be more 
and more localized around the origin as /io increases, as 
expected. 

Having checked that our solver performs quite well 
in all the above mentioned cases, it seemed tempting 
to further probe its capabilities with some other charge 
density distributions. In particular, we modelled a pla- 
nar capacitor (hence involving surface-like BC), and a 
cylindrical capacitor (wire-like BC). In the latter case, 
the input charge distribution was mimicked by a positive 
2D Gaussian distribution sharply peaked at the origin, 
together with a ring of negative Gaussian distributions 
centered at r = L/A. The relative amplitudes of the cen- 
tral and of the peripheral Gaussians were chosen to yield 
zero total charge. The solutions obtained are shown in 
Figs. 13-14, are definitely consistent with the intuitively 
expected behavior. 



VI. CONCLUSION 

We have presented a numerical method for the solution 
of Poisson's equation which can tackle any type of peri- 



Figure 11: Accuracy in the computation of the Hartree linear 
energy density - see Eq. (37) - in the case of wire-like boundary 
conditions (WBC) with monopolar charge density distribution 
(m is the order of the ISF, h the grid spacing). In our setup 
£ (exact) = 7 121H28333623614 hartree/bohr. 




Figure 12: Gaussian density charge distribution (red, dashed 
mesh) as a function of the isolated directions (x, y in our no- 
tation) and the corresponding electrostatic potential (black, 
solid mesh) evaluated for different values of the screening, 
namely /ig = {0, 0.01, 0.10, 1} bohr" 2 upon increasing bound- 
ary thickness. The charge density is implicitly periodic along 
z (wire-like BC). The amplitude of the plotted density charge 
distribution is multiplied by a factor 10 with respect to the 
actual value in order to improve the readability of the picture. 
Only half of the solution is drawn to highlight its profile. 



odicity, the presence of screening, non-orthorhombic ge- 
ometries and charged systems, showing that convergence 
to highly accurate results is attained with no adjustable 
parameter other than the (unavoidable) grid spacing. 

The charge density distribution and the electrostatic 
potential are both expressed in terms of plane waves 
along the periodic directions, and in terms of interpolat- 
ing scaling functions along the non-periodic directions. 
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Figure 13: Electrostatic potential generated by a pair of pla- 
nar charge distributions of opposite sign (positive in red/solid; 
negative in black/dashed) modelling a planar capacitor unlim- 
ited in the periodic directions. The piecewise planar behavior 
corresponds to the case with no screening (/io = 0), whereas 
the other solutions are obtained with /iq = {1, 10, 100} bohr -2 
upon increasing the boundary thickness. The potential is 
more and more localized around the capacitor's plates and 
falls rapidly to zero as fio is increased. Each curve is normalised 
to one to improve readability. 




Figure 14: Electrostatic potential generated by a cylindrical 
capacitor, periodic in the vertical direction. The different 
solutions correspond to /iq = {0, 1, 10, 100} bohr -2 upon in- 
creasing the boundary thickness. Each curve is normalised to 
one for sake of readability. Only half of the solution is drawn 
so as to highlight the potential profile. 



The latter representation proves to be very handy, in that 
the expansion coefficients of any continuous quantity are 
simply given by its values on a uniform grid. Moreover, 
rn-th order interpolating scaling functions preserve the 
matching between the moments built upon the discretized 
quantity with those of the continuous one. This occur- 
rence is particularly important in electrostatics, because 
of the interest in resolving the multipolar features of the 
electrostatic potential and other derived quantities. 
In our approach, the solution is obtained via the Green's 



function method. The (in principle) most expensive op- 
eration, namely the convolution of the Green's function 
with the input charge density distribution, becomes af- 
fordable by making use of highly-optimized 0(N log TV) 
FFT routines (zero- padded along the isolated directions), 
while the evaluation of the convolution kernel is carried 
out separately along the three spatial directions by ap- 
proximating the Green's function as a sum of Gaussian 
functions. 

Owing to the above mentioned advantages, our solver 
is suitable for intensive computer simulations of electronic 
structure and molecular dynamics, where the Poisson's 
equation has to be solved several times and it is important 
to limit the growth and propagation of numerical errors 
as much as possible, especially because other physical 
quantities are computed starting from the solution of the 
Poisson's equation. There are other contexts in computa- 
tional physics and chemistry in which relevant quantities 
are obtainable as convolutions involving the same Green's 
functions found in electrostatics. This is the case, for 
instance, of the exact exchange term within the general- 
izations of Kohn-Sham DFT employing orbital-dependent 
or hybrid functionals. Actually, our methodology features 
a level of generality which allows to address also problems 
well beyond electrostatics. 

As a possible outlook, we are working towards en- 
abling the computation of range-separated hybrid func- 
tionals, where the Coulomb potential is split into a 
long- and a short-range component. In the basic range- 
separated approach the Coulomb interaction is written as 
1/r = [erf(r/r ) + erfc(r/r )] /r, r being an adjustable 
length scale, although other options have been put for- 
ward (in which, for instance, the long- and short-range 
parts can be weighted differently). We are thus planning 
to model the range-separated Coulomb interaction within 
our framework. 

Our solver is already designed for taking full advantage 
of multi-core CPUs and GPUs, and is currently integrated 
in BigDFT, the sources of which are freely downloadable 
from http : / /inac . cea . f r/L_Sim/BigDFT/. The release 
of a stand-alone package is also envisaged for the near 
future. The details on the GPU acceleration and the per- 
formance of the solver in the context of massively parallel 
electronic structure computations will be described in a 
forthcoming paper. 
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